# MOC w/ linear tracking

from numpy import *
from scipy import *
from matplotlib.pyplot import * 
from operator import itemgetter
import time

class MOC(object):
    
    def __init__(self, nangle, delta, height, width, nrings, nsectors, npins, fuel_diameter):
        """
        All info pertains to one pin cell. npins recreates identical pin cells with one large grid.
        """
        self.nangle = nangle//2
        self.nsectors = nsectors
        self.nrings = nrings
        self.h = float(height)
        self.w = float(width)
        self.boundary = 'reflective'
        self.delta = delta
        self.regions = []
        self.rings = []
        self.fuel_diameter = fuel_diameter
        self.npins = npins
        self.pinCenters = []

        self.polarquad = Quadrature()
        self.polarquad.TabuchiYamamoto(3)

#############################################################################
#MAKE
#############################################################################
        
    def makeStuff(self):
        """Function to make tracks, segments, and solve for Keff"""
        tstart = time.time()
        self.fuel, self.water = Test()
        self.makeTracks()
        self.makeReflective()
        self.makePins()
        self.makeRegions()
        self.makeSegments()
        #self.solveKeff()
        self.plotSegments()
        tend = time.time()
        print 'Run time for ' + str(self.npins*(self.nrings+1)*self.nsectors) + ' regions: ' +  str(tend-tstart) + ' seconds'
        
    def makeSlopes(self):
        """Find the slopes for all tracks.
        nx - intersections on x-axis
        ny - intersections on y-axis
        nt - total intersections
        phi - angle
        m - slope
        n - segment length = n * (xEnd - xStart)
        """
        
        p = pi/self.nangle*(0.5 + arange(self.nangle))
        intv = vectorize(int)
        
        self.nx = intv(abs(self.w/self.delta*sin(p)) + 1)
        self.ny = intv(abs(self.h/self.delta*cos(p)) + 1)  
        self.nt = self.nx + self.ny                     
        self.phi = arctan(self.h/self.w*self.nx/self.ny)
        self.m = self.h/self.w*self.nx/self.ny
        self.n = sqrt(self.m**2 + 1)
        
        #correct angles with a negative slope
        for i, angle in enumerate(self.phi):
            if p[i] > pi/2:
                self.phi[i] = pi - self.phi[i]
                self.m[i] = -self.h/self.w*self.nx[i]/self.ny[i]
                self.n[i] = sqrt(self.m[i]**2 + 1)
       
    def makeAziWeights(self):
        """
        find the azimuthal weights
        """
        temp1 = concatenate((self.phi, [pi] + self.phi[0]))
        temp2 = concatenate((-self.phi[0:1], self.phi))
        x1 = 0.5*(temp1[1:]-temp1[:-1])
        x2 = 0.5*(temp2[1:]-temp2[:-1])
        self.wgta = (x1 + x2)/(2*pi)*self.deltas*2        
                    
    def makeDeltas(self):
        """
        xprime - separation along x-axis
        yprime - separation along y-axis
        deltas - track separation
        """
        self.xprime = self.w/self.nx
        self.yprime = self.h/self.ny
        self.deltas = self.xprime*sin(self.phi)
                    
    def makeTracks(self):
        #find the slopes, deltas, and weights for each angle
        self.makeSlopes()
        self.makeDeltas()
        self.makeAziWeights()
        
        #create lists of track slopes (m) and y-intercepts (b) for each angle
        self.tracks = []
        for i in range(self.nangle):
            b = zeros(self.nt[i])
            xin = zeros(self.nt[i])
            slope = zeros(self.nt[i])
            slope[::] = self.m[i]
            
            xin[:self.nx[i]] = self.xprime[i]*(0.5 + arange(self.nx[i]))
            b[:self.nx[i]] = -self.m[i]*xin[:self.nx[i]]
            
            if self.m[i] > 0:
                b[self.nx[i]:] = self.yprime[i]*(0.5 + arange(self.ny[i]))
            elif self.m[i] < 0:
                b[self.nx[i]:] = self.yprime[i]*self.nx[i] + self.yprime[i]*(0.5 + arange(self.ny[i]))
            
            # make a list of (m,b) and sort from min(b) to max(b)
            a = zip(slope, b)
            a = sorted(a, key = itemgetter(1))
                
            #create track2D objects for each track and store in list self.tracks   
            self.tracks.append([])
            for slope_cur,yint in a:
                newTrack = Track2D(slope_cur, yint, self.w, self.h, i, self.n[i], self.wgta[i])
                self.tracks[i].append(newTrack)
                
    def makeReflective(self):
        """ track_in         track #
            refl_in     is track flux going in the start (0) or end (1) of reflected track
            track_out        track #
            refl_out    is track flux going in the start (0) or end (1) of reflected track
       """
        
        for i in range(self.nangle//2):
            nx = self.nx[i]
            ny = self.ny[i]
            nt = self.nt[i]
            track_refl = self.tracks[-(i+1)]
            for j,track in enumerate(self.tracks[i]):
                if ny <= nx:
                    if j < ny:
                        #bottom to right
                        track.track_in                  = track_refl[nx - j - 1]
                        track.track_in.track_out        = track
                        track.refl_in                   = 1
                        track.track_in.refl_out         = 0
                       
                        track.track_out                 = track_refl[nx + j]
                        track.track_out.track_out       = track
                        track.refl_out                  = 1
                        track.track_out.refl_out        = 1
                    
                    elif j < nx:
                        #bottom to top
                        track.track_in                  = track_refl[nx - j - 1]
                        track.track_in.track_out        = track
                        track.refl_in                   = 1
                        track.track_in.refl_out         = 0
                                                
                        track.track_out                 = track_refl[nt+ny-j-1]
                        track.track_out.track_in        = track
                        track.refl_out                  = 0
                        track.track_out.refl_in         = 1
                        
                    else:
                        #left to top
                        track.track_in                  = track_refl[j - nx]
                        track.track_in.track_in         = track
                        track.refl_in                   = 0
                        track.track_in.refl_in          = 0
                        
                        track.track_out                 = track_refl[nt+ny-j-1]
                        track.track_out.track_in        = track
                        track.refl_out                  = 0
                        track.track_out.refl_in         = 1
                        
                        
                else:
                    if j < nx:
                        #bottom to left
                        track.track_in                  = track_refl[nx-j-1]
                        track.track_in.track_out        = track
                        track.refl_in                   = 1
                        track.track_in.refl_out         = 0

                        track.track_out                 = track_refl[nx+j]
                        track.track_out.track_out       = track
                        track.refl_out                  = 1
                        track.track_out.refl_out        = 1
                        
                        
                    elif j < ny:
                        # left to right
                        track.track_in                  = track_refl[j - nx]
                        track.track_in.track_in         = track
                        track.refl_in                   = 0
                        track.track_in.refl_in          = 0
                        
                        track.track_out                 = track_refl[nx+j]
                        track.track_out.track_out       = track
                        track.refl_out                  = 1
                        track.track_out.refl_out        = 1
                        
                    else:
                        # left to top
                        track.track_in                  = track_refl[j - nx]
                        track.track_in.track_in         = track
                        track.refl_in                   = 0
                        track.track_in.refl_in          = 0

                        track.track_out                 = track_refl[nt+ny-j-1]
                        track.track_out.track_in        = track
                        track.refl_out                  = 0
                        track.track_out.refl_in         = 1
                        
    def makeRegions(self):
        """
        Create Region2D objects for each region within the pin cell. The fuel pin is comprised of rings with the outermost ring being Ring 1. Ring 0 is tecnically the water.  Sectors correspond to angular pie slices.  The water sectors are numbered first (regions 0-nsectors) followed by fuel pin going from outermost ring towards the center.
        """
        npins = self.npins
        nrings = self.nrings
        nsectors = self.nsectors
        
        #total volume of the fuel pin and each ring within the fuel pin
        Vfuel = pi*(float(self.fuel_diameter)/2)**2
        Vring = Vfuel/(nrings)
        for i in range(nrings):
            # find ring radius
            ring = nrings - i
            if i == 0:
                radius = sqrt(Vring/pi)
                self.rings.append(radius)
            else:
                radius = sqrt(Vring/pi + self.regions[0].radius**2)
                self.rings.append(radius)
                #create Region2D object for each sector of the fuel rings.
            for pin in self.pinCenters:
                for sector in range(nsectors):
                    newRegion = Region2D(radius, ring, nrings, sector, nsectors, self.fuel, pin)
                    newRegion.material_name = 'fuel'
                    newRegion.volume = Vring/nsectors
                    self.regions.insert(0,newRegion)
            
        #acreate Region2D object for each sector of the water
        for pin in self.pinCenters:            
            for sector in range(nsectors):
                newRegion = Region2D(None, 0, nrings, sector , nsectors, self.water, pin)
                newRegion.material_name = 'water'
                newRegion.volume = self.findWaterVolume(newRegion)
                self.regions.insert(0, newRegion)
            
         #sort regions based on region number   
        self.regions = sorted(self.regions, key = lambda Region2D: Region2D.region)
        
    def makeSegments(self):
        """First, all interesections of a track with rings and sectors are found.  Segments are then created by stitching together intersection points and stored with their corresponding track.  A pointer is included within the Segment2D object to link it to a region
        """
        #loop over tracks
        for track in self:
            m = track.m
            b = track.b
            
            newSegments = []
                
            #make segments on the cell boundaries
            trackIn = track.findIn()
            trackOut = track.findOut()
            newSegments.append(trackIn)
            newSegments.append(trackOut)
            
            # find all ring intersections
            ringSegments = self.findRingIntersections(track)
            # find all sector intersections
            sectorSegments = self.findSectorIntersections(track)
            # find all pin cell intersections
            cellSegments = self.findCellIntersections(track)
            
            #add ring and sector segments to newSegments list
            for k in ringSegments:
                newSegments.append(k)
            for k in sectorSegments:
                newSegments.append(k)
            for k in cellSegments:
                newSegments.append(k)
                
            #sort segments from xmin to xmax as we are tracking from left to right    
            sortedSegments = sorted(newSegments)
                
            #remove duplicates
            sortSegments = []
            SegmentOld = -2
            for v in sortedSegments:
                if abs(v - SegmentOld) > 10e-10:
                    sortSegments.append(v)
                SegmentOld = v
            
            #loop over segments in track   
            for j in range(len(sortSegments)-1):
                start = sortSegments[j]
                end = sortSegments[j+1]
                middle = float(start + end)/2.0000001
                # use the middle point on the track to find what region it is in
                region = self.findRegion(track, middle)
                ystart = m*start + b
                yend = m*end + b
                #stitch together start and endpoint
                newSegment = Segment2D(region,start,end,ystart,yend)
                newSegment.length = track.n *(end - start)
                track.segments.append(newSegment)
                self.regions[region].segments.append(newSegment)
                self.regions[region].volume_calc += newSegment.length*track.weight
                
        #correct region lengths to account for different track/volume ratios for each region         
        for region in self.regions:
            region.factor = region.volume/region.volume_calc
        for track in self:
            for segment in track:
                segment.length *= self.regions[segment.region].factor
                
                
    def makePins(self):
        """
        make pin centers.  There only only be 1,4,9,16,25... pins.
        """
        npins = self.npins
        npins_side = int(round(sqrt(npins)))
        if npins == 1:
            xBounds = [0, self.w]
            yBounds = [0, self.h]
            newPin = Pin2D(self.w/2,self.h/2,0, xBounds, yBounds)
            self.pinCenters.append(newPin)
        else:
            for i in range(npins_side):
                for j in range(npins_side):
                    pin = int(round(sqrt(npins)*i + j))
                    xBounds = [i*self.h/npins_side, (i+1)*self.w/npins_side]
                    yBounds = [j*self.h/npins_side, (j+1)*self.h/npins_side]
                    xCenter = (i+.5)*self.w/npins_side
                    yCenter = (j+.5)*self.h/npins_side
                    newPin = Pin2D(xCenter, yCenter, pin, xBounds, yBounds)
                    self.pinCenters.append(newPin)
            
        

###########################################################################
#SOLVE
###########################################################################
                
    def solve(self):
        """Solve the backwards characteristic form of the transport
        equation over the domain with a fixed source (assumed to be
        flat in each region). We start by assuming that the fluxes on
        all the boundaries are zero and then iterate until they
        converge."""

        # Precompute exponential factor and set incoming fluxes to zero

        for track in self:
            track.flux_in = zeros((2,3))
            for s in track:
                sigma = self.regions[s.region].sigma_t
                sintheta = self.polarquad.sin
                s.factor = 1.0 - exp(-sigma*s.length/sintheta)                    
        

        # Iterate on incoming fluxes on domain boundary
        max_iterations = 100
        for i_boundary in range(max_iterations):
            #print("Inner iteration {0}".format(i_boundary))

            # Initialize volume flux to zero
            for r in self.regions:
                r.flux = 0

            # Loop over polar angles
            for j, (sintheta, wgtp) in enumerate(self.polarquad):

                # Loop over each track
                for track in self:
                    # Need to divide total weight for each track by 2
                    # since we are tracking forwards and backwards
                    weight = 4*pi*track.weight*wgtp*sintheta

                    # Loop over each segment in track
                    flux = track.flux_in[0,j]
                    if i_boundary == 0:
                        flux = 0
                    for s in track:
                        r = self.regions[s.region]
                        delta = (flux - r.source/r.sigma_t)*s.factor[j]
                        r.flux += delta*weight
                        flux -= delta

                    # Transfer flux to outgoing track
                    track.track_out.flux_in[track.refl_out, j] = flux

                    # Now traverse this track in reverse for reflective
                    if self.boundary == "reflective":
                        flux = track.flux_in[1,j]
                        if i_boundary == 0:
                            flux = 0
                        for s in track.segments[::-1]:
                            r = self.regions[s.region]
                            delta = (flux - r.source/r.sigma_t)*s.factor[j]
                            r.flux += delta*weight
                            flux -= delta

                        # Transfer flux to incoming track
                        track.track_in.flux_in[track.refl_in, j] = flux


            # Add in source and normalize to volume
            for r in self.regions:
                r.flux /= 2
                r.flux = 4*pi*r.source/r.sigma_t + r.flux/(r.sigma_t*r.volume)

            # Save source in fuel for Dancoff factor
            if i_boundary == 0:
                for region in self.regions:
                    if region.material_name == 'fuel':
                        region.flux0 = region.flux

            # Check if converged
            fluxv = array([r.flux for r in self.regions])
            if i_boundary > 0:
                if all(abs((fluxv - fluxvold)/fluxvold) < 1e-10):
                    return
            fluxvold = fluxv
                    
                    
    def solveKeff(self, max_outer = 1000):
        """Using the fixed source kernel from above, iterate on the
        scattering and fission sources until they converge (outer
        iterations)."""

        # Initial guess
        S_old = {}
        S = {}
        for r in self.regions:
            S_old[r.region] = 1
        k_old = 1.0

        # Loop over outer iterations
        for i_outer in range(max_outer):
            self.counter = i_outer
            # Inner iteration
            for r in self.regions:
                scatter = r.sigma_s*r.flux/(4*pi)
                r.source = 1/(4*pi*k_old)*S_old[r.region] + scatter
            self.solve()

            # Outer iteration
            for r in self.regions:
                S[r.region] = r.nusigma_f*r.flux
            keff = k_old*sum(S.values())/sum(S_old.values())
            
            print("i = {0}, keff = {1}".format(i_outer, keff))
            if abs(keff - k_old) < 1e-7:
                break

            # Update source
            k_old = keff
            for r in self.regions:
                S_old[r.region] = S[r.region]
 
#########################################################################
#FIND
#########################################################################
                
    def findWaterVolume(self, region):
        """
        find the volume of each sector of the water.  Since the water region is not symmetric, we cannot create a general equation for the volume of each slice.  The code is only able to calculate the volume for 1,2,4,8,or 16 sector divisions
        """
        cell_side = self.h
        if self.npins != 1:
            cell_side = self.h/(self.npins//2)
        vol_total = cell_side**2 - pi*(self.fuel_diameter/2)**2
        if self.nsectors == 1 or 2 or 4 or 8:
            return vol_total/self.nsectors
        elif self.nsectors == 16:
            if region == 1 or 2 or 5 or 6 or 9 or 10 or 13 or 14:
                return self.w**2/8*(1-tan(pi/8))-pi/16*(fuel_diameter/2)**2
            else:
                return self.w**2/8*tan(pi/8)-pi/16*(fuel_diameter/2)**2
            
    def findRingIntersections(self, track):
        """Use the quadratic equation to find the intersections of the track with the rings. The general form of a circle (x - xCenter)**2 + (y - yCenter)**2 = radius**2 is used to represent each ring. 
        """ 
        
        h = self.h
        w = self.w
        
        newRingSegments = []
        for pin in self.pinCenters:
            for k in range(len(self.rings)):
                #create quadratic parameters
                x = pin.x
                y = pin.y
                a = 1 + track.m**2
                b = 2*track.m*track.b - 2*track.m*y - 2*x
                c = x**2 + y**2 + track.b**2 - 2*track.b*y - self.rings[k]**2
                c = c.real
                
                quad = b**2 - 4*a*c
                
                #find and store interesections
                if quad == 0:
                    int1 = -b/(2*a)
                    newRingSegments.append(int1)
                elif quad > 0:
                    int1 = (-b - sqrt(quad))/(2*a)
                    int2 = (-b + sqrt(quad))/(2*a)
                    newRingSegments.append(int1)
                    newRingSegments.append(int2)
                
        return newRingSegments
        
    def findSectorIntersections(self, track):
        """find intersections of a track with the sector divisions. Each sector division passes through the center of the pin cell and is represented with a slope and y-intercept. 
        """
        sectorSegments = []
        if self.nsectors != 1:
            #loop over sectors
            for pin in self.pinCenters:
                for i in range(self.nsectors):
                    m = self.regions[i].m1
                    b = pin.y - m*pin.x
                    if m != track.m:
                        #find intersection point
                        xint = (b - track.b)/(track.m - m)
                        yint = track.m * xint + track.b
                        #find out if intersection point is within the geometry
                        if xint > pin.xBCs[0] and xint < pin.xBCs[1] and yint > pin.yBCs[0] and yint < pin.yBCs[1]:
                            sectorSegments.append(xint)
           
        return sectorSegments
    
    def findCellIntersections(self, track):
        """
        find the intersections of a track and the cell boundaries
        """
        CellSegments = []
        
        for pin in self.pinCenters:
            yint_1 = track.m*pin.xBCs[0] + track.b
            yint_2 = track.m*pin.xBCs[1] + track.b
            xint_1 = (pin.yBCs[0] - track.b)/track.m
            xint_2 = (pin.yBCs[1] - track.b)/track.m
            
            if yint_1 >= pin.yBCs[0] and yint_1 <= pin.yBCs[1]:
                CellSegments.append(pin.xBCs[0])
            if yint_2 >= pin.yBCs[0] and yint_2 <= pin.yBCs[1]:
                CellSegments.append(pin.xBCs[1])
            if xint_1 >= pin.xBCs[0] and xint_1 <= pin.xBCs[1]:
                CellSegments.append(xint_1)
            if xint_2 >= pin.xBCs[0] and xint_2 <= pin.xBCs[1]:
                CellSegments.append(xint_2)
                
        return CellSegments
        
        
    
    def findRegion(self, track, x):
        """
        Given an x value for a point on a track, find out what region that point is in.  The procedure is based on first finding what ring the point is in followed by what sector divisions it lies between.  With this info, the point (and therefore the segment the point is on) is assigned a region.
        """
        m = track.m
        b = track.b
        h = self.h
        w = self.w 
        
        y = m*(x) + b
        
        #find out what pin cell the point is in
        radius = 10e5
        for pin_test in self.pinCenters:
            radius_temp = sqrt(abs(y-pin_test.y)**2 + abs(x-pin_test.x)**2)
            if radius_temp < radius:
                pin = pin_test
                radius = radius_temp
        
        #find what ring the point is in
        ring = 0
        for i in range(len(self.rings)):
            c1 = len(self.rings) - i - 1
            if radius < self.rings[c1]:
                ring = self.nrings - c1
        
        #find angle between radius line and y = height/2
        if y >= pin.y:
            ptangle = arccos((x - pin.x)/radius)
        else:
            ptangle = pi + arccos((pin.x-x)/radius)
            
        #find what sector the point is in
        for i in range(self.nsectors):
            if ptangle < self.regions[i].angle2 and ptangle > self.regions[i].angle1:
                sector = i
        
        #calculate region number based on pin, sector, and ring number
        region = pin.number*(self.nsectors*(self.nrings+1)) + self.nsectors*ring + sector
        return region

#########################################################################
#PLOT
#########################################################################
    
    def plotSegments(self):
        """
        Assign all segments of a region the same color.  Plot all segments.
        """
        k_x = []
        k_y = []
        b_x = []
        b_y = []
        c_x = []
        c_y = []
        r_x = []
        r_y = []
        y_x = []
        y_y = []
    

        for track in self:
            for segment in track:
                if segment.region % 5 == 0:
                    k_x.append(segment.segment[0])
                    k_y.append(segment.segment[1])
                elif segment.region % 5 == 1:
                    b_x.append(segment.segment[0])
                    b_y.append(segment.segment[1])
                elif segment.region % 5 == 2:
                    c_x.append(segment.segment[0])
                    c_y.append(segment.segment[1])
                elif segment.region % 5 == 3:
                    r_x.append(segment.segment[0])
                    r_y.append(segment.segment[1])
                elif segment.region % 5 == 4:
                    y_x.append(segment.segment[0])
                    y_y.append(segment.segment[1])

        t0 = time.time()
        xlist = []
        ylist = []
        for xends,yends in zip(k_x,k_y):
            xlist.extend(xends)
            xlist.append(None)
            ylist.extend(yends)
            ylist.append(None)
        plot(xlist,ylist, 'k-')
        
        xlist = []
        ylist = []
        for xends,yends in zip(b_x,b_y):
            xlist.extend(xends)
            xlist.append(None)
            ylist.extend(yends)
            ylist.append(None)
        plot(xlist,ylist, 'b-')
        
        xlist = []
        ylist = []
        for xends,yends in zip(c_x,c_y):
            xlist.extend(xends)
            xlist.append(None)
            ylist.extend(yends)
            ylist.append(None)
        plot(xlist,ylist, 'c-')
        
        xlist = []
        ylist = []
        for xends,yends in zip(r_x,r_y):
            xlist.extend(xends)
            xlist.append(None)
            ylist.extend(yends)
            ylist.append(None)
        plot(xlist,ylist, 'r-')
        
        xlist = []
        ylist = []
        for xends,yends in zip(y_x,y_y):
            xlist.extend(xends)
            xlist.append(None)
            ylist.extend(yends)
            ylist.append(None)
        plot(xlist,ylist, 'y-')
        
        t1 = time.time()
        segments = 0
        for track in self:
            segments += len(track.segments)
        print 'Plotting time for ' + str(segments) + ' segments: ' +  str(t1-t0) + ' seconds'
        
        
    def plotReflective(self, angle, trackNum):
        """
        Plot a track and the tracks it intersects with at the boundaries.
        """
        track = self.tracks[angle][trackNum]
        Track2D.plot(track)
        Track2D.plot(track.track_in)
        Track2D.plot(track.track_out)
        
    def plotTracks(self):
        """
        plot all tracks
        """
        for i in range(self.nangle):
            for j in self.tracks[i]:
                Track2D.plot(j)   
    
    def plotRegion(self, region):
        """
        plot all segments for a given region
        """
        for track in self:
            for segment in track:
                if segment.region == region:
                    plot(segment.segment[0],segment.segment[1], 'k-')                
  
######################################################################
#MISCELLANEOUS 
######################################################################
                    
    def printData(self):
        for track in self:
                for segment in track:
                    print track.angle, segment.region, segment.segment
    
    def __iter__(self):
        for i in range(self.nangle):
            for track in self.tracks[i]:
                yield track
                
                        
class Segment2D(object):
    """
    Every segment is represented as a Segment2D object with attributes specifying the region it is in and its start and end points.
    """
    
    def __init__(self, region, xStart, xEnd, yStart, yEnd):
        self.region = region
        self.segment = [[xStart, xEnd],[yStart, yEnd]]

class Region2D(object):
    """
    Every region is represented as a Region2D object.
    angle1 - angle (clockwise from horizontal) of minimum angle
    angle2 - angle (clockwise from horizontal) of maximum angle
    m1 - slope of minimum angle sector boundary
    b1 - y-intercept of mimimum angle sector boundary
    m2 - slope of maximum angle sector boundary
    b2 - y-intercept of maximum angle sector boundary
    volume - geometric volume of region
    volume_calc - volume calculated from segment lengths
    sigma_t, sigma_s, nusigma_f - total, scattering and fission xs's
    segments - segments within region
    flux - flux in region
    source - source in region
    radius - outer radius boundary of region
    region - region number
    """
    
    def __init__(self, radius, ring, nrings, sector, nsectors, material, pin):
        self.radius = radius
        self.sigma_t = material.sigma_t
        self.sigma_s = material.sigma_s
        self.nusigma_f = material.nusigma_f
        self.volume_calc = 0.
        self.volume = 0.
        self.flux = 0.
        self.source = material.source
        self.segments = []
        self.pin = pin
        self.pinNumber = pin.number

        x = pin.x
        y = pin.y
        self.region = pin.number*(nsectors*(nrings+1)) + nsectors*ring + sector
        self.angle1 = 2*pi*sector/nsectors
        if abs(self.angle1 - pi/2) < 10e-10:
            self.m1 = 10e10
            self.b1 = y - self.m1*x
        elif abs(self.angle1 - 3*pi/2) < 10e-10:
            self.m1 = 10e10
            self.b1 = y - self.m1*x
        else:
            self.m1 = tan(self.angle1)
            self.b1 = y - self.m1*x
        if abs(self.m1) < 10e-10:
            self.m1 = 0
            self.b1 = y
        self.angle2 = 2*pi*(sector+1)/nsectors
        if abs(self.angle2 - pi/2) < 10e-10:
            self.m2 = 10e10
            self.b2 = y - self.m2*x
        elif abs(self.angle2 - 3*pi/2) < 10e-10:
            self.m2 = 10e10
            self.b2 = y - self.m2*x
        else:
            self.m2 = tan(self.angle2)
            self.b2 = y - self.m2*x
        if abs(self.m2) < 10e-10:
            self.m2 = 0
            self.b2 = y
            
            
class Track2D(object):
    """
    Object used to represent each track. All tracks are represented as a line of the form y = m*x + b. 
    h - geometry height
    w - geometry width
    m - slope
    b - y intercept
    n - segment length = n * (xEnd - xStart)
    wgta - track weight
    angle - angle number (not actual angle)
    track_in - pointer to incoming track
    track_out - pointer to outgoing track
    segments - list of segments within track. sorted from xmin to xmax
    """
    
    def __init__(self, m, b, w, h, angle = None, n = None, wgta = None):
        self.h        = h
        self.w        = w
        self.m        = m
        self.b        = b
        self.n        = n
        self.weight   = wgta
        self.angle    = angle
        self.track_in = 0
        self.track_out= 0
        self.segments = []     
            
    def plot(self, linespec = 'k-'):
        x1 = self.findIn()
        y1 = self.m*x1 + self.b
        x2 = self.findOut()
        y2 = self.m*x2 + self.b
        plot([x1,x2],[y1,y2],linespec)
        

    def findIn(self):
        """
        finds where track enters geometry
        """
        if self.m > 0:
            yout = self.m*self.w + self.b
            if yout < self.h:
                if self.b < 0:
                    return -self.b/self.m
                elif self.b > 0:
                    return 0
            elif yout > self.h:
                if self.b < 0:
                    return -self.b/self.m
                elif self.b > 0:
                    return 0
        
        elif self.m < 0:
            yout = self.m*self.w + self.b
            if yout > 0:
                if self.b > self.h:
                    return (self.h-self.b)/self.m
                elif self.b < self.h:
                    return 0
            elif yout < 0:
                if self.b > self.h:
                        return (self.h-self.b)/self.m
                if self.b < self.h:
                    return 0
                
    def findOut(self):
        """
        finds where track exits geometry
        """
        if self.m > 0:
            yout = self.m*self.w + self.b
            if yout < self.h:
                if self.b < 0:
                    return self.w
                elif self.b > 0:
                    return self.w
            elif yout > self.h:
                if self.b < 0:
                    return (self.h-self.b)/self.m
                elif self.b > 0:
                    return(self.h-self.b)/self.m
                    
        elif self.m < 0:
            yout = self.m*self.w + self.b
            if yout > 0:
                if self.b > self.h:
                    return self.w
                elif self.b < self.h:
                    return self.w
            elif yout < 0:
                if self.b > self.h:
                        return -self.b/self.m
                if self.b < self.h:
                    return -self.b/self.m   
                
    def __iter__(self):
        for segment in self.segments:
            yield segment
            
class Pin2D(object):
    
    def __init__(self, x, y, pin, xBounds, yBounds):
        self.x = x
        self.y = y
        self.number = pin
        self.xBCs = xBounds
        self.yBCs = yBounds
        self.source = 0.
        self.flux = 0.
        self.flux0 = 0.
    
                
class Quadrature(object):
    """
    Polar quadrature used for method of characteristics. Right now the
    only options are to use the Leonard optimized quadrature of the
    Tabuchi-Yamamoto quadrature.

    Attributes:
      sin    = sin of the polar angle (measured from z-axis)
      weight = corresponding polar weight

    Methods:
      TabuchiYamamoto = set quadrature to a TY quadrature 
      Leonard         = set quadrature to Leonard optimized quadrature
      __init__        = initialize quadrature
      __iter__        = enable ability to iterate over pairs of (sin, weight)
    """

    def __init__(self):
        self.sin = []
        self.weight = []

    def __iter__(self):
        for i in range(len(self.sin)):
            yield (self.sin[i], self.weight[i])

    def __len__(self):
        return len(self.sin)

    def TabuchiYamamoto(self, n):
        """
        Set the polar angles and weights according to the
        Tabuchi-Yamamoto polar quadrature derived in A. Yamamoto et
        al., "Derivation of Optimum Polar Angle Quadrature Set for the
        Method of Characteristics Based on Approximation Error for the
        Bickley Function," J. Nucl. Sci. Tech., 44 (2), pp. 129-136
        (2007).
        """

        if n == 1:
            self.sin = array([0.798184])
            self.weight = array([1.0])
        elif n == 2:
            self.sin = array([0.363900, 0.899900])
            self.weight = array([0.212854, 0.787146])
        elif n == 3:
            self.sin = array([0.166648, 0.537707, 0.932954])
            self.weight = array([0.046233, 0.283619, 0.670148])

    def Leonard(self, n):
        """
        Set the polar angles and weights according to the Leonord
        Optimized polar quadrature from A. Leonard and C. McDaniel,
        "Optimal Polar Angles and Weight for the Characteristics
        Method," Trans. Am. Nucl. Soc., 73, pp. 172-173 (1995) as well
        as R. Le Tellier and A. Hebert, "Anisotropy and Particle
        Conservation for Trajectory-Based Deterministic Methods,"
        Nucl. Sci. Eng., 158, pp. 28-39 (2008).
        """

        if n == 2:
            self.sin = array([0.273658, 0.865714])
            self.weight = array([0.139473, 0.860527])
        elif n == 3:
            self.sin = array([0.099812, 0.395534, 0.891439])
            self.weight = array([0.017620, 0.188561, 0.793819])


class Material(object):
    """
    object that defines the properties of a material
    """
    def __init__(self, name, sigma_t, sigma_s, nusigma_f, source = 0.):
        self.name = name
        self.sigma_t = sigma_t
        self.sigma_s = sigma_s
        self.nusigma_f = nusigma_f
        self.source = source

######################################################################
#PROBLEMS
######################################################################
        
def Test():
    # Adjust fuel cross-sections
    fuel = Material('fuel', 4.52648699E-01,3.83259177E-01, 9.94076580E-02)

    # Adjust moderator cross-sections
    water = Material('water', 8.41545641E-01, 8.37794542E-01, 0)
    
    return fuel, water

def Problem2():

    # Adjust fuel cross-sections
    fuel = Material('fuel', 4.52648699E-01,3.83259177E-01, 9.94076580E-02)

    # Adjust moderator cross-sections
    water = Material('water', 8.41545641E-01, 8.37794542E-01, 0)
    
    nrings = 1
    nsectors = 1
    npins = 4
    pins_side = 2
    tstart = time.time()
    q = MOC(16,.05,1.27*pins_side,1.27*pins_side,nrings,nsectors,npins, .8)
    q.makeTracks()
    q.makeReflective()
    q.fuel, q.water = fuel, water
    q.makePins()
    q.makeRegions()
    q.makeSegments()
    q.plotSegments()
    q.solveKeff() 
    tend = time.time()
    print 'Run time: ' +  str(tend-tstart) + ' seconds'
    
def Problem1Part1():
    """Create table of fuel/coolant flux ratio as a function of ray
    spacing and number of azimuthal angles for fuel xs = 1.0 and
    coolant xs = 0.1
    """
    tstart = time.time()
    fuel = Material('fuel', 1, 0, 0, 1./(4*pi))
    
    water = Material('water', 0.1, 0, 0, 0.)
    
    rings = 1
    sectors = 1
    pins = 4
    pins_side = 2
    for nangle in [16,32,64,128]:#,256,512]:
        for spacing in [0.10, 0.05, 0.02]:#, 0.01]:
            q = MOC(nangle, spacing, 1.27*pins_side, 1.27*pins_side, rings, sectors, pins, .8)
            q.makeTracks()
            q.makeReflective()
            q.fuel, q.water = fuel, water 
            q.makePins()
            q.makeRegions()
            q.makeSegments()
            q.solve()
            for pin in q.pinCenters:
                water_flux = 0.
                for region in q.regions:
                    if region.material_name == 'water' and region.pinNumber == pin.number:
                        water_flux += region.flux
                    elif region.material_name == 'fuel' and region.pinNumber == pin.number:
                        print(pin.number, region.region, nangle, spacing, region.flux/water_flux)
                    
    tend = time.time()
    print 'Run time: ' +  str(tend-tstart) + ' seconds'
                
            
def Problem1Part2():
    """Make six tables of Dancoff factor (fuel xs = infinity) and
    coolant xs = 0, 1, and infinity as a function of ray spacing and
    number of azimuthal angles
    """
    tstart = time.time()
    fuel = Material('fuel', 1e5, 0, 0, 1./(4*pi))
    
    nrings = 1
    nsectors = 1
    npins = 9
    npins_side = 3
    for sigma_t in [1e-5, 1, 1e5]:
        water = Material('water', sigma_t, 0, 0 , 0.)
        for nangle in [16,32,64,128]:#,256,512]:
            for spacing in [0.10, 0.05, 0.02]:#, 0.01]:
                q = MOC(nangle, spacing, 1.27*npins_side, 1.27*npins_side, nrings, nsectors, npins, .8)
                q.makeTracks()
                q.makeReflective()
                q.fuel, q.water = fuel, water 
                q.makePins()
                q.makeRegions()
                q.makeSegments()
                q.solve()
                for pin in q.pinCenters:
                    for region in q.regions:
                        if region.pinNumber == pin.number and region.material_name != 'water':
                            if region.material_name == 'fuel':
                                pin.source += region.source
                                pin.flux += region.flux
                                pin.flux0 += region.flux0
                                pin.sigma_t = region.sigma_t
                    dancoff = (1 - (pin.source/pin.sigma_t*4*pi - pin.flux)/(pin.source/pin.sigma_t*4*pi - pin.flux0))
                    print(pin.number, sigma_t, nangle, spacing, dancoff)
                        
    tend = time.time()
    print 'Run time: ' +  str(tend-tstart) + ' seconds'
                